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ABSTRACT.— Predicting available habitat and the fine scale distribution of an endemic butterfly lizard, 
Leiolepis ocellata, is important for conservation since its populations are experiencing threats from anthropogenic 
activities. We conducted field surveys in northern Thailand and detected L. ocellata at 15 of 16 surveyed 
localities. We then used Maxent to predict its potential distribution based on presence-only data collected from 
field surveys and museum specimens and predictor variables (19 bioclimatic variables). The most optimal model 
confirmed its geographic distribution in northern Thailand but absent from other regions in the country supported 
by high AUC (0.875+0.076). The most important predictor variable was isothermality followed by temperature 
annual range and precipitation seasonality. Its occurrence was corresponded with known sites for the species. One 
further site, from which an old specimen was collected and recorded in 1967, was not predicted by the model and 
disappeared during field surveys. Most predicted distribution range located outside protected areas, which may 
potentially disturb L. ocellata populations. Our results presented updated information on current occurrence and 
predicted distribution of L. ocellata in northern Thailand, which is useful for developing conservation strategy of 
this species. 
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INTRODUCTION 


The eyed butterfly lizard Leiolepis 
ocellata Peters, 1971 is an endemic species 
in northern Thailand. This species was first 
described as ZŁ. belliana ocellata, a 
subspecies of L. belliana (Hardwicke et 
Gray, 1827) by Peters (1971), and later it 
has been followed by subsequent authors 
(Tiedemann et al., 1994; Manthey, 2010; 
Chan-ard et al., 2015). But it is currently 
elevated as a full species based on the 


difference in number of chromosome 
(Arunyavalai, 2003), distinct external 
morphology (Arunyavalai, 2003), and 


distinguishable phylogenetic relationship 


(Pauwels and Chimsunchart, 2007; Das, 
2010; Grismer et al., 2014). 

Leiolepis ocellata ecologically serves as 
an intermediate prey for important higher- 
level predators such as snakes and raptors. It 
becomes food source for rural people in 
northern Thailand and likely overhunted in 
some areas. Other human activities such as 
wildfire for agricultural purposes and urban 
expansion could also deteriorate habitats of 
L. ocellata (Arunyavalai, 2003). Without 
any protection act in Thailand and also 
listed as “Least Concern” by the 
International Union Conservation of Nature 
(IUCN) Red List, L. ocellata is potentially 
decreasing and may lead to local extinct. 
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Baseline information especially on 
distribution of this endemic taxon is an 
urgent priority for conservation planning. 
Northern Thailand is topographically 
mountainous with many large mountain 
ranges separating lowlands. These 
intermontane lowlands are important 
habitats for L. ocellata. It becomes a 
challenge to evaluate the distribution of L. 
ocellata since an information on its 
distribution is scarce and not up to date. 
Largest source of information came from 
Arunyavalai (2003) who reported L. 
ocellata in Phayao, Nan, Lampang, and Tak 
Province. We opportunistically took 
photographs of L. ocellata in protected areas 
and forest park of the Forest Industry 
Organization of Thailand (FIO) such as Doi 
Suthep-Pui National Park in Chiang Mai 
Province, Doi Wiang La Wildlife Sanctuary 
in Mae Hong Son Province, and Thung 
Kwian and Mae Chang Forest Park in 
Lampang Province. Reassessment of its 
distribution may reveal new unexplored 
populations that will be a novel resource for 
future research and resource management. 
Species distribution modeling becomes 
useful for creating distribution map based 


on known localities and background 
environmental conditions (Guisan and 
Zimmermann, 2000). The models that 


require both presence and absence data on 
species of interest could be developed with 
variety of statistical methods (Phillips et al., 
2006). However, true absence data is 
sometimes unavailable for the scattered data 
sources such as old reports and museum 
specimens that usually described only 
localities of capture (Phillips et al., 2006; 
Ponder et al., 2001). Maximum entropy 
modeling (Maxent) which requires only 
presence data, has been developed to cope 
with this limitation (Phillips et al., 2006). 
Maxent estimates habitat suitability by 


finding the most uniform distribution 
(maximum entropy) across the areas of 
interest given environmental constraints. 
Climatic variables, in particular 
temperature, are considered major 
environmental constraints for ectothermic 
animals (Buckley et al., 2008; Shine, 2005; 
Spellerberg, 1971) and influenced 
distribution modeling of reptiles (e.g. Gül, 
2015; Javed et al., 2017). Thus, fluctuations 
in temperature and perhaps other 
environmental climatic variables across 
diverse topography in northern Thailand are 
expected to play a significant role in 
distribution modeling of L. ocellata. 

The objective of this study was to survey 
and predict the distribution of L. ocellata 
based on presence-only data obtained from 
field surveys and museum collection and 
bioclimatic variables. We limited our 
analysis to Thailand where we could follow 
the previous records and perform field 
surveys. We expected that the model could 
predict distribution range of this endemic 
species which is crucial for conservation 
and spatial planning. 


MATERIALS AND METHODS 


Field surveys 

We selected localities for our survey 
localities based on historical records 
including Arunyavalai (2003), reports from 
the Forest Industry Organization of Thailand 
(FIO) and the Department of National Parks, 
Wildlife, and Plant Conservation (DNP), 
museum records (four specimens deposited 
at Thailand Natural History Museum; 
voucher No. THNHM15653, THNHM16101, 
THNHM 23215, and THNHM 25974), local 
reports and interviews with local people. 
Surveys were conducted at 16 known 
localities (Table 1) using visual encounter 
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TABLE 1. Localities data of L. ocellata based on field surveys and museum records 

















Detection 
Localities Provinces Evidences during field 
surveys 
Ban Tak District Tak Local report (2017) Yes 
Mae Wa National Park, Mae Phrik District* Tak Arunyavalai (2003) Yes 
Ban Mae La Mao, Mae Sod District Tak THNHM15653 (1967) No 
Si Satchanalai Forest Park, Si Satchanala District Sukhothai FIO report (2013) and Not surveyed 
THNHM23215 (2013) 
Doi Pha Klong National Park, Long District* Phrae Local report (2017) Yes 
Pha Sing Subdistrict, Mueang Nan District Nan Local report (2018) Yes 
Wiang Sa District Nan Arunyavalai (2003) Yes 
Mae Chai Subdistrict , Mae Chai District Phayao Local report (2017) Yes 
Mae Peum National Park, Pa Deat District* Chiang Rai Local report (2017) Yes 
Mae Chang Forest Park, Mae Mo District Lampang FIO report (2015) and Not surveyed 
THNHM25974 (2015) 
Mueang Lampang District Lampang Local report (2017) Yes 
Chae Son Subdistrict, Mueang Pan District Lampang Local report (2017) Yes 
Thung Kwian Forest Park, Hang Chat District Lampang FIO report (2013) Yes 
Mae Li Forest Park, Li District Lamphun FIO report (2016) Yes 
Ban Khilek, Mae Rim District Chiang Mai THNHM16101 (1970) Yes 
Doi Suthep-Pui National Park, Maeung Chiang Mai District* Chiang Mai Local report (2017) Yes 
Mae Chaem Forest Park, Mae Chaem District Chiang Mai FIO report (2015) Yes 
Doi Wiang La Wildlife Sanctuary, Khun Yuam District* Mae Hong Son Local report (2018) Yes 
* = protected area 
survey technique (VES). Most surveys were 2017). This resolution approximately 


situated in forest parks and national parks 
with some localities located on private 
property and school grounds. The observers 
walked along existing trails in those 
localities in order to minimize disturbance 
and carefully searched for butterfly lizards. 
The surveys were carried out during 9:00 to 
15:00 h on sunny days during the dry season 
from October 2017 to May 2018. To 
increase the probability of detection, we 
conducted 3-5 surveys for each locality on 
separated days. Coordinates of present 
localities were recorded. Two points of 
detection located < 1 km in distance apart 
were considered the same point. 
Environmental variables 

All bioclimatic variables (Supplementary 
-Table S1) were downloaded from WorldClim 
database (www.worldclim.org) at a resolution 
of 30 sec (~1 km2) under monthly average 
conditions for 1970-2000 (Fick and Hijmans, 


corresponded with the estimated home range 
of Leiolepis (approximately 800 m2; 
Arunyavalai, 2003). In Thailand, the initial 
prediction of Maxent showed that 
distribution of L. ocellata did not appeared 
in other regions but northern Thailand. We 
thus presented only the map of northern 
Thailand. The climate layers were clipped to 
the north of Thailand using QGIS. All layers 
were then converted into ASCH files. 
Species distribution modeling 

Presence data and environmental layers 
were used to generate probability maps 
predicting the potential distribution of L. 
ocellata in Maxent v.3.4.1 (Phillips et al., 
2006). Presence data was obtained from our 
field surveys and specimen data (four 
localities with authenticated museum 
specimens). To select the most suitable 
predictor variables, 19 bioclimatic variables 
were initially considered in the first 
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modeling using default setting of the 
Maxent (regularization multiplier = 1, max 
number of background points = 10000, and 
1 replicate with cross-validate replicated run 
type). The selection was based on percent 
contribution (>10%), permutation importance 
(>10%), and Pearson correlation (coefficient 
Ir| < 0.80; e.g. Khanum et al., 2013). After 
the initial model was evaluated, only three 
variables, isothermality (bio3), temperature 
annual range (bio7), and precipitation 
seasonality (bio15), were extracted and used 
for further analyses. Isothermality (bio3), 
which determines day-night temperature 
oscillation relative to annual temperature 
oscillations (summer-winter, O’Donnell and 
Ignizio, 2012), was calculated to percent. A 
value of 100% indicated equality of day- 
night temperature range and annual 
temperature range, while values < 100% 


indicated that day-night temperature 
oscillations were smaller than annual 
temperature oscillations. | Temperature 


annual range (bio7) was calculated by 
subtracting the minimum temperature of the 
coldest month from the maximum 
temperature of the warmest month 
(O’Donnell and Ignizio, 2012). Precipitation 
seasonality (bio15) is a measure of variation 
in monthly precipitation totals over the 
course of the year (O’Donnell and Ignizio, 


2012). 

Twenty models were generated using 
different combinations of these three 
variables and different regularization 


multiplier values (1-3). These experiments 
were run under 10 replications of cross- 
validation replicated run type, the output 
format as “Cloglog”, 10 percentile training 
presence as threshold rule, and other 
parameters at its default setting. Model 
performance was evaluated from areas 
under the curve (AUC) of a receiver 
operating characteristic (ROC) plot. The 


experiment with the highest average test 
AUC was considered as the most optimal 
model. The contribution of each 
environmental variable was then assessed by 
two methods: 1) percent contribution and 
permutation importance and 2) the jackknife 


method. Response curves were also 
generated to observe how each 
environmental variable affects Maxent 


prediction. The most optimal model was 
projected with presence data of L. ocellata. 
Further, the raster results were polygonized 
to calculate a total area and the area located 
within protected areas. Projection and 
polygonization were generated using QGIS. 


RESULTS 


Leiolepis ocellata were highly active in 
the mornings of sunny days but less active 
and usually remained in their burrows on 
rainy days. Habitat types of the surveyed 
sites consisted of bare grounds, agricultural 
areas, and dry dipterocarp forests. The areas 
were encompassed by mountain ranges and 
covered with loose soil, across an elevation 
range of 140 to 646 meters above sea level. 
We detected L. ocellata at 35 points in 15 
discrete localities (Table 1). Of these, there 
were five localities located in protected 
areas. We also detected additional L. 
ocellata populations at four sites that had 
never been documented but were informed 
by local people, which were Mae Peum 
National Park (Chiang Rai Province), Mae 
Chai District (Phayao Province), Doi Pha 
Klong National Park (Phrae Province) and 
Khun Yuam District (Mae Hong Son 
Province). The only locality where we failed 
to detect L. ocellata was in Mae Sod 
District, Tak Province. Nevertheless, we 
included this locality based on a museum 
specimen (THNHM15653, collected in 
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TABLE 2. Model selection summary of L. ocellata. (Grey bar represents the most optimal model.) 











Combination of variables Number of variables Regularization multiplier value Average test AUC 
bio3, bio7, biol5 3 1 0.872+0.095 
bio3, bio7, biol5 3 1.5 0.875+0.085 
bio3, bio7, biol5 3 2 0.875+0.081 
bio3, bio7, biol5 3 2.5 0.875+0.081 
bio3, bio7, biol5 3 3 0.875+0.076 
bio3, bio7 2 1 0.854+0.077 
bio3, bio7 2 1.5 0.857+0.073 
bio3, bio7 2 2 0.858+0.070 
bio3, bio7 2 2.5 0.857+0.068 
bio3, bio7 2 3 0.855+0.067 
bio3, biol5 2 1 0.844+0.062 
bio3, biol 5 2 1.5 0.848+0.060 
bio3, biol5 2 2 0.851+0.062 
bio3, biol5 2 2.5 0.854+0.066 
bio3, biol5 2 3 0.855+0.068 
bio7, bio1l5 2 1 0.869+0.106 
bio7, bio1l5 2 1.5 0.871+0.104 
bio7, bio1l5 2 2 0.871+0.101 
bio7, bio1l5 2 2.5 0.872+0.098 
bio7, bio1l5 2 3 0.862+0.092 














1967) as a confirmed site for our model 
prediction. As aresult, a total of 39 points of 
L. ocellata— 35 points from field surveys 
and four points from museum records— 
were used in model prediction. 

Optimization of Maxent generated 20 
distribution models supported by average 
test AUC values (Table 2). The model with 
a combination of three variables (bio3, bio7, 
biol5) showed the highest AUC (0.875) 
when regularization multiplier values were 
1.5-3. Thus, the most optimal model with 
lowest standard deviation (0.875+0.076) 
was used as a representative model for 
predicting the potential distribution of L. 
ocellata (Table 2; Fig1). 

The areas with high probability of 
occurrence were predicted in northern 
Thailand, especially in the middle of the 
region (Fig. 1), and absent from other 
regions in the country as predicted by initial 
modeling. The distribution modeling mostly 
corresponded with surveyed localities and 





museum records with the sole exception of 
the record from Mae Sod District, Tak 
Province (Fig. 1). This locality specified by 
one museum specimen was neither predicted 
by this model nor detected during field 
surveys. Percent mutation, permutation 
importance, and jackknife procedure 
showed that isothermality (bio3) was the 
most important predictor variable (Table 3; 
Fig. 2) followed by temperature annual 
range (bio7) and precipitation seasonality 
(biol5). The probability of occurrence 
decreased with increasing of isothermality 
and precipitation seasonality (Fig. 3A, C) 
but increased with increasing temperature 
annual range (Fig. 3B). 

After polygonization, the most optimal 
model estimated 10,710 km2 as a total 
potential distribution area of L. ocellata 
suitable (Fig. 4). Of this, 2,879 km2 located 
across 32 protected areas in northern 
Thailand (ca. 27% of total predicted area). 
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FIGURE 1. Predicted distribution model based on three selected predictor variables (bio3, 7, and 15) and 
presence information based on surveys (©) and museum specimens (starsyr) of L. ocellata. Mae Sod District, 
Tak Province, was a locality obtained from museum specimen, but not found in predicted map and survey (0O). 


DISCUSSION 


Improved and more detailed information 
on distribution of poorly known or endemic 
species is fundamental for conservation 
(Burgman et al., 2005; Gül, 2015; Javed et 
al., 2017). Field surveys and museum 
records enabled us to generate the first 
species distribution modeling for L. 
ocellata. The use of presence-only data 
should be interpreted with caution since an 
imperfect prediction may result from limited 
algorithms of the model (Elith and Graham, 
2009). Nevertheless, we reduced this 


imperfection by optimization to find the 
most reliable model with the average AUC 
> 0.8 (Eskildsen et al., 2013; Manel et al., 
2001; Swets, 1988). The most optimal 
model (AUC = 0.875+0.076) revealed the 
areas with high probability of occurrence of 
L. ocellata in northern Thailand and absent 
from other region in the country. Of the 
three predictor variables selected for 
performing the model, isothermality (bio3) 
had the strongest influence on L. ocellata 
distribution. Isothermality indicating 
fluctuation of monthly temperature ranges 
throughout the year was negatively 
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TABLE 3. Estimated relative contributions of each environmental variable. The values are average over 10 











replications. 
Variable Percent contribution Permutation importance 
Isothermality (bio3) 64.4 61.1 
Temperature annual range (bi07) 24.2 15.1 








Precipitation seasonality (bio15) 11.4 23.9 
2 (A) Jackknife of regularized training gain 
e bio15 Without variable 
> With only variable ™ 
£ bio3 u With all variables I 
E i E 
5 bio7 
Z 
w 
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Without variable 
With only variable ™ 
J777 With all variables M 


0.80 0.85 


FIGURE 2. Jackknife test of the final model presenting environmental variable contributions to (A) 
regularized training gain and (B) AUC of L. ocellata. The values are average over 10 replications. 


correlated with the probability of occurrence 
of L. ocellata. This suggested L. ocellata 
may prefer habitats with fluctuated 
temperature throughout the year as found in 
northern Thailand. 

Within northern Thailand, low probability 
of occurrence was mostly found in the areas 
along mountain ranges in which the temperature 
is low throughout the year (occasionally 
dropping below 2 °C; www.tmd.go.th). As 
environmental temperature is vital for 
ectotherms (Angilletta et al., 2002; Lutter- 
schmidt and Hutchison, 1997), this 
temperature may be lower than the critical 
thermal minimum of lizards (approximately 
2.2 - 9.8 °C; Spellerberg, 1971) and could 


limit distribution of L. ocellata. Moreover, 
the highest mountain range in Thailand, 
Thanon Thong Chai Range, isolates the 
lowland in Mae Hong Son Province in the 
extreme westernmost part from other areas. 
Mae Hong Son Province is a poorly 
explored area in which biological 
information on butterfly lizards and others 
reptile species is still lacking. Our first 
report and model prediction of L. ocellata in 
this region suggested the possibility of many 
more unexplored reptiles that may exist and 
require comprehensive investigation. 

The model prediction well corresponded 
with historical records and field surveys 
except for one locality obtained from a 
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FIGURE 3. Response curves presenting relationships between selected environmental variables (x axis) and 
probability of presence of L. ocellata (y axis). (A) isothermality (bio3), (B) temperature annual range (bi07), 
and (C) precipitation seasonality (biol5). Each curve represents a different Maxent model created using only 


the corresponding variable. 


museum specimen (THNHM15653) 
collected in Mae Sod District in 1967. Mae 
Sod District was not suitable for L. ocellata 
as predicted by the model and confirmed 
from field surveys. We carefully conducted 
field surveys and did not find the butterfly 
lizards nor their burrows in this and adjacent 
areas. Because important information such 
as name of the collector was missing, it was 
possible that mislabeling of a document 
might have occurred for Mae Sod District. 
Baseline knowledge on species 
distribution in Thailand is imperative but 
usually unorganized. In this study, we 
gathered several sources of information for 
predicting potential distribution of L. 
ocellata. The model created a useful insight 


on the species biogeography and 
conservation. We found that predicted range 
of L. ocellata mostly located outside 
protected areas such as in private 
agricultural fields and urban zones, which 
are highly modified. Their populations in 
those areas are, thus, potentially threatened 
by land modification. For instance, in 
Mueang Chiang Mai District, Chiang Mai 
Province, we found L. ocellata only in a 
small part of Doi Suthep foothill despite its 
predicted suitable habitat covered most part 
of the district. Urban expansion in Chiang 
Mai Province may relegate its habitats 
preventing population establishment. 
Furthermore, the protected areas containing 
predicted range were mostly dry deciduous 
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FIGURE 4. Predicted distribution model covered by protected areas in northern Thailand. 


forests and currently have been severely 
devastated by wildfire. We strongly 
suggested species distribution modeling 
could be applied for conservational purposes 
of L. ocellata including the establishment of 
additional protected areas and prevention of 
those areas from urbanization and wildfire. 
We also recommend that organized and 
updated species distribution data could be 
useful for conservation planning of others 
threatened/endemic species, which will 
require consistent and careful monitoring. 
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Supplementary material 


SUPPLEMENTARY TABLE S1. Bioclimatic variables used in modeling 














Name Code 
Annual mean temperature biol 
Mean diurnal range (mean of monthly (max temp — min temp)) bio2 
Isothermality (bio2/bio7) (*100) bio3 
Temperature seasonality (standard deviation * 100) bio4 
Max temperature of warmest month bios 
Min temperature of coldest month bio6 
Temperature annual range (bio5-bio6) bio7 
Mean temperature of wettest quarter bio8 
Mean temperature of driest quarter bio9 
Mean temperature of warmest quarter biolO 
Mean temperature of coldest quarter bioll 
Annual precipitation biol2 
Precipitation of wettest month biol3 
Precipitation of driest month biol4 
Precipitation seasonality (coefficient of variation) biol5 
Precipitation of wettest quarter biol6 
Precipitation of driest quarter biol7 
Precipitation of warmest quarter biol8 
Precipitation of coldest quarter biol9 











